clear
tic

%spec determines which specification to use.  1 uses CRRA 1 specification of
%Stock and Wright (2000), while 2 adds moments based on the interest rate,
%to repredoce the results reported in the supplementary appendix of 
%Andrews (2017).

spec=2;

%Set nominal coverage for confidence sets: nominal coverage is 1-alpha
alpha=0.05;
%alpha=0.1;

%Sets the values of gamma_min to be used in two-step confidence set
%construction
gamma_min_vec=ones(3,1)*.05;
%gamma_min_vec=ones(3,1)*.5;

load crra;
n=length(w);
if spec==1
Z=[ones(n-2,1) w(2:n-1)-mean(w(2:n-1)) r(2:n-1)-mean(r(2:n-1))];
data=[w(3:n) r(3:n) Z];
elseif spec==2
Z=[ones(n-2,1) w(2:n-1)-mean(w(2:n-1)) r(2:n-1)-mean(r(2:n-1)) ...
    ones(n-2,1) w(2:n-1)-mean(w(2:n-1)) rf(2:n-1)-mean(rf(2:n-1))];
data=[w(3:n) r(3:n) rf(3:n) Z];    
end

%Set grids of parameter values to consider
 delta_grid=[0:0.0025: 1.1];
 eta_grid=[-6:0.05:200];
% delta_grid=[0:0.025: 1];
% eta_grid=[-6:0.5:150];

theta_grid_array={delta_grid,eta_grid};

%Set which functions of parameter to consider when constructing confidence
%sets
F1=eye(2);
F2=[1 0; 0 0];
F3=[0 0; 0 1];
Farray=F1;
Farray(:,:,2)=F2;
Farray(:,:,3)=F3;

%Set which moment conditions and weighting matricies to use
if spec==1
giveMoments='HS_moment_condition_generalized_spec1';
giveWeight='Efficient_weight_spec1';
elseif spec==2
    giveMoments='HS_moment_condition_generalized_spec2';
    giveWeight='Efficient_weight_spec2';
end
%giveWeight='Homosked_weight';
giveVarianceEstim='NW_hac';
[CS_R_stack, CS_N_stack, gamma_hat_vec]=Twostep_CS_calc(alpha,gamma_min_vec,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim);

%Report results: for the confidence set on the full (two-dimensional)
%parameter vector reports a scatter plot, while for
for n=1:size(Farray,3)
    if n==1
        b1=delta_grid;
        b2=eta_grid;
        edges={b1,b2};
        
        points=CS_R_stack(n);
        points=cell2mat(points);
        
        figure(n);
        
        N=hist3(points,'Edges',edges);
        CS_R=N>0;
        colormap([1 1 1;.7 .7 .7])
        
        points=CS_N_stack(n);
        points=cell2mat(points);
        N=hist3(points,'Edges',edges);
        CS_N=N>0;
        colormap([1 1 1;.7 .7 .7])

        figure(1)
        subplot(1,2,1)
        contourf(b2,b1,CS_R,1,'k','Linewidth',2)
        xlabel('\eta','Fontsize',20)
        ylabel('\delta','Fontsize',20)
        title('CS_R','Fontsize',12)
        subplot(1,2,2)
        contourf(b2,b1,CS_N,1,'k','Linewidth',2)
        xlabel('\eta','Fontsize',20)
        ylabel('\delta','Fontsize',20)
        title('CS_N','Fontsize',12)
    elseif n==2
        points=CS_R_stack(n);
        points=cell2mat(points);
        points=unique(points);
        CS_R_inclusion=ismember(delta_grid,points);
        
        points=CS_N_stack(n);
        points=cell2mat(points);
        points=unique(points);
        CS_N_inclusion=ismember(delta_grid,points);
        
        figure(n);
        subplot(1,2,1)
        plot(delta_grid,CS_R_inclusion);
        subplot(1,2,2)
        plot(delta_grid,CS_N_inclusion);
    else
        points=CS_R_stack(n);
        points=cell2mat(points);
        points=unique(points);
        CS_R_inclusion=ismember(eta_grid,points);
        
        points=CS_N_stack(n);
        points=cell2mat(points);
        points=unique(points);
        CS_N_inclusion=ismember(eta_grid,points);
        
        figure(n);
        subplot(1,2,1)
        plot(eta_grid,CS_R_inclusion);
        subplot(1,2,2)
        plot(eta_grid,CS_N_inclusion);
    end
end